clear all

t=0.5; num=50; k=2;dt=t/num;
time=linspace(0,t,num+1);
exact=zeros(size(time)); numerical=exact; 
exactfun=@(s) 2/t*t^2/(k-1)^2/pi^2*sin((k-1)*pi/t*s);
vals=exactfun(time);
for ii=2:num+1;
    exact(ii)=quadl(exactfun,0,time(ii));
    numerical(ii)=midptrule(vals(1:ii),dt);
end
hold on
plot(time,abs(exact-numerical))
%plot(time,numerical,'r')
